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Inhomogeneous temporal processes, like those appearing in human communications, neuron spike 
trains, and seismic signals, consist of high-activity bursty intervals alternating with long low-activity 
periods. In recent studies such bursty behavior has been characterized by a fat-tailed inter-event time 
distribution, while temporal correlations were measured by the autocorrelation function. However, 
these characteristic functions are not capable to fully characterize temporally correlated heteroge- 
nous behavior. Here we show that the distribution of the number of events in a bursty period serves 
as a good indicator of the dependencies, leading to the universal observation of power-law distribu- 
tion in a broad class of phenomena. We find that the correlations in these quite different systems 
can be commonly interpreted by memory effects and described by a simple phenomenological model, 
which displays temporal behavior qualitatively similar to that in real systems. 



In nature there are various phenomena, from earth- 
quakes [1 to sunspots [2 and neuronal activity, that 
show temporally inhomogeneous sequence of events, in 
which the overall dynamics is determined by aggregate 
effects of competing processes. This happens also in hu- 
man dynamics as a result of individual decision mak- 
ing and of various kinds of correlations with one's so- 
cial environment. These systems can be characterized 
by intermittent switching between periods of low ac- 
tivity and high activity bursts [4-6 , which can appear 
as a collective phenomena similar to processes seen in 
self-organized criticality [7HT6]. In contrast with such 
self-organized patterns intermittent switching can be de- 
tected at the individual level as well (see Figjl]), as seen 
for single neuron firings and for earthquakes at a single 
location, both showing inhomogeneous temporal patterns 

[9H11111I1IIH1E2]. 

Further examples of bursty behavior at the individual 
level have been observed in the digital records of human 
communication activities through different channels [H 
[23H26] . Over the last few years different explanations 
have been proposed about the origin of inhomogeneous 
human dynamics [H [22l [27], including the single event 
level [28^, and about the impact of circadian and weekly 
fluctuations [7]. Moreover, by using novel technology of 
Radio Frequency ID's, heterogeneous temporal behavior 
was observed in the dynamics of face-to-face interactions 
[3Ql [3Tj . This was explained by a reinforcement dynamics 
[32l |33] driving the decision making process at the single 
entity level. 

For systems with discrete event dynamics it is usual to 
characterize the observed temporal inhomogeneities by 
the inter-event time distributions, P(tie), where tie — 
tij^i — ti denotes the time between consecutive events. 
A broad P{tie) O |23[ |34] reflects large variability in 
the inter-event times and denotes heterogeneous tempo- 
ral behavior. Note that P{tie) alone tells nothing about 
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the presence of correlations, usually characterized by the 
autocorrelation function, A(r), or by the power spectrum 
density. However, for temporally heterogeneous signals 
of independent events with fat-tailed P(tie) the Hurst ex- 
ponent can assign false positive correlations [l] together 
with the autocorrelation function (see Supplementary In- 
formation) . To understand the mechanisms behind these 
phenomena, it is important to know whether there are 
true correlations in these systems. Hence for systems 
showing fat-tailed inter-event time distributions, there is 
a need to develop new measures that are sensitive to cor- 
relations but insensitive to fat tails. 




20d:0h 2h 4h 6h 8h lOh 12h 14h 16h 18h 20h 22h 21d:0h 



FIG. 1. Activity of single entities with color-coded inter-event 
times, (a): Sequence of earthquakes at a single location (b): 
Firing sequence of a single neuron (c): Outgoing mobile phone 
call sequence of an individual. Shorter the time between the 
consecutive events darker the color. 



In this paper we define a new measure that is capable 
of detecting whether temporal correlations are present, 
even in the case of heterogeneous signals. By analyzing 
the empirical datasets of human communication, earth- 
quake activity, and neuron spike trains, we observe uni- 
versal features induced by temporal correlations. In the 
analysis we establish a close relationship between the ob- 
served correlations and memory effects and propose a 
phenomenological model that implements memory driven 
correlated behavior. 
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I. RESULTS 



A. Correlated events 



A sequence of discrete temporal events can be inter- 
preted as a time-dependent point process, where 
= 1 at each time step ti when an event takes place, 
otherwise X{ti) =0. To detect bursty clusters in this bi- 
nary event sequence we have to identify those events we 
consider correlated. The smallest temporal scale at which 
correlations can emerge in the dynamics is between con- 
secutive events. If only X{t) is known, we can assume two 
consecutive actions at ti and t^+i to be related if they fol- 
low each other within a short time interval, t^+i —ti<At 
[28l [36] . For events with the duration di this condition 
is slightly modified: t^+i — {ti -\- di) < At. 

This definition allows us to detect bursty periods, de- 
fined as a sequence of events where each event follows the 
previous one within a time interval At. By counting the 
number of events, that belong to the same bursty pe- 
riod, we can calculate their distribution P{E) in a signal. 
For a sequence of independent events, P{E) is uniquely 
determined by the inter-event time distribution P{tie) as 
follows: 



P{E = n) = ^ P{t^e)dt,ej ^ l P{t^e)dt,ej 

(1) 

for n > 0. If the measured time window is finite, the 
integral J^^ P{tie)dtie < 1 and P{E = n) - a-^^-^), 

otherwise if P{tie)dtie = 1 then P{E = n) = 
(for related numerical results see SI). Consequently for 
any finite independent event sequence the P{E) distri- 
bution decays exponentially even if the inter-event time 
distribution is fat-tailed. Deviations from this exponen- 
tial behavior indicate correlations in the timing of the 
consecutive events. 



1. Bursty sequences in human communication 

To check the scaling behavior of P{E) in real systems 
we focused on outgoing events of individuals in three se- 
lected datasets: (a) A mobile-call dataset from a Euro- 
pean operator; (b) Text message records from the same 
dataset; (c) Email communication sequences [24 (for de- 
tailed data description see Materials and Methods). For 
each of these event sequences the distribution of inter- 
event times measured between outgoing events are shown 
in Figj2] (left bottom panels) and the estimated power- 
law exponent values are summarized in Table |lj To ex- 
plore the scaling behavior of the autocorrelation func- 
tion, we took the averages over 1,000 randomly selected 
users with maximum time lag of r = 10^. In Figj2]a 
and b (right bottom panels) for mobile communication 
sequences strong temporal correlation can be observed 
(for exponents see Table |l|. The power- law behavior in 



A{r) appears after a short period denoting the reaction 
time through the corresponding channel and lasts up to 
12 hours, capturing the natural rhythm of human activi- 
ties. For emails in Figj2jc (right bottom panels) long term 
correlation are detected up to 8 hours, which refiects a 
typical office hour rhythm (note that the dataset includes 
internal email communication of a university staff). 
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TABLE I. Characteristic exponents of the (ol) autocorrelation 
function, (/3) bursty number, (7) inter-event time distribution 
functions and v memory functions (see SI) calculated in dif- 
ferent datasets and for the model study. 

The broad shape of P{tie) and A{r) functions con- 
firm that human communication dynamics is inhomoge- 
neous and displays non-trivial correlations up to finite 
time scales. However, after destroying event-event cor- 
relations by shuffling inter-event times in the sequences 
(see Materials and Methods) the autocorrelation func- 
tions still show slow power-law like decay (empty symbols 
on bottom right panels), indicating spurious unexpected 
dependencies. This clearly demonstrates the disability of 
A{r) to characterize correlations for heterogeneous sig- 
nals (for further results see SI). However, a more effec- 
tive measure of such correlations is provided by P{E). 
Calculating this distribution for various At windows, we 
find that the P{E) shows the following scale invariant 
behavior 



P{E) - E- 



(2) 



for each of the event sequences as depicted in the main 
panels of Fig j2] Consequently P{E) captures strong tem- 
poral correlations in the empirical sequences and it is re- 
markably different from P{E) calculated for independent 
events, which, as predicted by ([T]), show exponential de- 
cay (empty symbols on the main panels). 

Exponential behavior of P{E) was also expected from 
results published in the literature assuming human com- 
munication behavior to be uncorrelated [271 EH] • How- 
ever, the observed scaling behavior of P{E) offers direct 
evidence of correlations in human dynamics, which can 
be responsible for the heterogeneous temporal behavior. 
These correlations induce long bursty trains in the event 
sequence rather than short bursts of independent events. 

We have found that the scaling of the P{E) distribu- 
tion is quite robust against changes in At for an extended 
regime of time- window sizes (Figj2|. In addition, the 
measurements performed on the mobile-call sequences in- 
dicate that the P{E) distribution remains fat-tailed also 
when it is calculated for users grouped by their activity. 
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FIG. 2. The characteristic functions of human communication event sequences. The P{E) distributions with various At 
time-window sizes (main panels), P{tie) distributions (left bottom panels) and average autocorrelation functions (right bottom 
panels) calculated for different communication datasets. (a) Mobile-call dataset: the scale- invariant behavior was characterized 
by power-law functions with exponent values a c^i 0.5, /3 ^ 4.1 and 7 c^i 0.7 (b) Almost the same exponents were estimated 
for short message sequences taking values a ^ 0.6, /3 ^ 3.9 and 7 ^ 0.7. (c) Email event sequence with estimated exponents 
a :^ 0.75, /3 2.5 and 7 = 1.0. Empty symbols assign the corresponding calculation results on independent sequences. Lanes 
labeled with s, m, h and d are denoting seconds, minutes, hours and days respectively. 



Moreover, the observed scaling behavior of the character- 
istic functions remains similar if we remove daily fluctu- 
ations (for results see SI). These analyses together show 
that the detected correlated behavior is not an artifact of 
the averaging method nor can be attributed to variations 
in activity levels or circadian fluctuations. 



2. Bursty periods in natural phenomena 

As discussed above, temporal inhomogeneities are 
present in the dynamics of several natural phenomena, 
e.g. in recurrent seismic activities at the same loca- 
tion pHTTj (for details see Materials and Methods and 
SI). The broad distribution of inter-earthquake times in 
Figjsja (right top panel) demonstrates the temporal in- 
homogeneities and the long tail of the autocorrelation 
function (right bottom panel) assigning long-range tem- 
poral correlations. Counting the number of earthquakes 
belonging to the same bursty period with At = 2... 32 
hours window sizes, we obtain a broad P{E) distribution 
(see Figjsja main panel), as observed earlier in commu- 
nication sequences, but with a different exponent value 
(see in Table [l|). Note that the presence of long bursty 
trains in earthquake sequences were already assigned to 
long temporal correlations by measurements using con- 
ditional probabilities [38l [39] . 

Another example of naturally occurring bursty behav- 
ior is provided by the firing patterns of single neurons 
(see Materials and Methods). The recorded neural spike 
sequences display correlated and strongly inhomogeneous 
temporal bursty behavior, as shown in Figjsjb. The dis- 
tributions of the length of neural spike trains are found 



to be fat-tailed and indicate the presence of correlations 
between consecutive bursty spikes of the same neuron. 



3. Memory process 

In each studied system (communication of individu- 
als, earthquakes at a given location, or single neurons) 
qualitatively similar behavior was detected as the sin- 
gle entities performed independent events or they passed 
through longer correlated bursty cascades. These cas- 
cades can be viewed as the result of building up some 
stress in the system (even at the single component level) 
e.g. the accumulation of important tasks in human com- 
munication, mechanical stress preceding an earthquake 
or integrated stimuli to a neuron. Our hypothesis is that 
each bursty period is related to a single excitation that 
triggers a cascade relieving this stress, which would nat- 
urally explain the correlation between the events within 
the same train. This "one cascade-one task" hypothe- 
sis can be relevant for human communication, where we 
speculate that a single task execution can occasionally 
trigger a series of events responsible for the appearance 
of long bursty trains. This assumption is supported by 
the observation that bursty nodes communicate predom- 
inantly with only one of their neighbors 23 SO] , indicat- 
ing that a bursty period maybe linked to a single task. 
For earthquakes we can make the same conjecture, since 
a bursty period of earthquakes at a given time and loca- 
tion are most likely related to the same seismic activity. 
In case of neurons, the firings take place in bursty spike 
trains when the neuron receives excitatory input and its 
membrane potential exceeds a given potential threshold 
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FIG. 3. The characteristic functions of event sequences of 
natural phenomena. The P(E) distributions of correlated 
event numbers with various At time- window sizes (main pan- 
els), P{tie) distributions (right top panels) and average auto- 
correlation functions (right bottom panels), (a) One station 
records of Japanese earthquake sequences from 1986 to 1998. 
The functional behavior is characterized by the fitted power- 
law functions with corresponding exponents a :^ 0.3, /3 2.5 
and 7 ^ 0.7. Inter-event times for P{tie) were counted with 10 
second resolution, (b) Firing sequences of single neurons with 
2 millisecond resolution. The corresponding exponents take 
values as a ^ 1.0, /3 :^ 2.3 and 7 1.1. Empty symbols as- 
sign the calculation results on independent sequences. Lanes 
labeled with ms, s, m, h, d and w are denoting milliseconds, 
seconds, minutes, hours, days and weeks respectively. 



[4T] . The spikes fired in a single train are correlated since 
they are the result of the same excitation and their firing 
frequency is coding the amplitude of the incoming stimuli 

m- 

The correlations taking place between consecutive 
bursty events can be interpreted as a memory process, 
allowing us to calculate the probability that the entity 
will perform one more event within a At time frame af- 
ter it executed n events previously in the actual cascade. 
This probability can be written as: 



p{n) 



E 



P{E) 



(3) 



Therefore the memory function, p(n), gives a different 
representation of the distribution P{E). The p{n) calcu- 
lated for the mobile call sequence are shown in Figj4ja for 
trains detected with different window sizes. Note that in 



empirical sequences for trains with size smaller than the 
longest train, it is possible to have p{n) = 1 since the 
corresponding probability would be P{E = n) = 0. At 
the same time due to the finite size of the data sequence 
the length of the longest bursty train is limited such that 
p{n) shows a finite cutoff. 

We can use the memory function to simulate a se- 
quence of correlated events. If the simulated sequence 
satisfies the scaling condition in ([2| we can derive the 
corresponding memory function by substituting ^ into 
(|3|, leading to: 



p{n) 



with the scaling relation (see SI): 
f3 = u^l. 



(4) 



(5) 



In order to check whether ([5| holds for real systems 
and whether the memory function in Q describes cor- 
rectly the memory in real processes we compare it to 
a memory function extracted from an empirical P{E) 
distributions. We selected the P{E) distribution of the 
mobile call dataset with At = 600 second and derived 
the corresponding p{n) function. The complement of the 
memory function, 1 —p(n)^ is presented in FigJIJb where 
we show the original function with strong finite size ef- 
fects (grey dots) and the same function after logarithmic 
binning (black dots). 

Taking equation Q we fit the theoretical memory 
function to the log-binned empirical results using least- 
squares method with only one free parameter, v. We 
find that the best fit offers an excellent agreement with 
the empirical data (see Fig 4jb and also FigJIJa) with 
V = 2.971±0.072. This would indicate p 3.971 through 
([5|, close to the approximated value /3 4.1, obtained 
from directly fitting the empirical P{E) distributions in 
the main panel of Figj2]a (for fits of other datasets see SI). 
In order to validate whether our approximation is cor- 
rect we take the theoretical memory function p{n) of the 
form Q with parameter u = 2.971 and generate bursty 
trains of 10^ events. As shown in Figjsjc, the scaling of 
the P{E) distribution obtained for the simulated event 
trains is similar to the empirical function, demonstrating 
the validity of the chosen analytical form for the memory 
function. 



Model study 

As the systems we analysed are of quite different na- 
ture, from physics (earthquakes) to social (human com- 
munication) and biological (neuron spikes) systems, find- 
ing a single mechanistic model to fit them all is virtually 
impossible. Therefore, our goal is to define a phenomeno- 
logical model that captures common features of the ob- 
served dynamics and see how these features are related 
to each other. 
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FIG. 4. Empirical and fitted memory functions of tfie mobile call sequence (a) Memory function calculated from the mobile 
call sequence using different At time windows, (b) 1 —p{n) complement of the memory function measured from the mobile call 
sequence with At = 600 second and fitted with the analytical curve defined in equation Q with u = 2.971. Grey symbols are 
the original points, while black symbols denotes the same function after logarithmic binning, (c) P{E) distributions measured 
in real and in modeled event sequences. 



4' Reinforcement dynamics with memory 

We assume that the investigated systems can be de- 
scribed with a two-state model, where an entity can be 
in a normal state A, executing independent events with 
longer inter-event times, or in an excited state 5, per- 
forming correlated events with higher frequency, corre- 
sponding to the observed bursts. To induce the inter- 
event times between the consecutive events we apply a 
reinforcement process based on the assumption that the 
longer the system waits after an event, the larger the 
probability that it will keep waiting. Such dynamics 
shows strongly heterogeneous temporal features as dis- 
cussed in [32j |33]. For our two-state model system we 
define a process, where the generation of the actual inter- 
event time depends on the current state of the system. 
The inter-event times are induced by the reinforcement 
functions that give the probability to wait one time unit 
longer after the system has waited already time tie since 
the last event. These functions are defined as 



fA,B{tie) = 



tie 



(6) 



where jua and jub control the reinforcement dynamics in 
state A and 5, respectively. These functions follow the 
same form as the previously defined memory function in 
Q and satisfy the corresponding scaling relation in ([5|. 
If jiiA <^ Mb the characteristic inter-event times at state 
A and B become fairly different, which induces further 
temporal inhomogeneities in the dynamics. The actual 
state of the system is determined by transition proba- 
bilities shown in FigjSjb, where to introduce correlations 
between consecutive excited events performed in state B 
we utilize the memory function defined in equation Q. 

To be specific, the model is defined as follows: first the 
system performs an event in a randomly chosen initial 
state. If the last event was in the normal state A, it waits 
for a time induced by fA{tie)^ after which it switches to 
excited state B with probability tt and performs an event 
in the excited state, or with probability 1 — tt stays in the 
normal state A and executes a new normal event. In the 



excited state the inter-event time for the actual event 
comes from fsitie) after which the system decides to 
execute one more excited event in state B with a prob- 
ability p{n) that depends on the number n of excited 
events since the last event in normal state. Otherwise it 
switches back to a normal state with probability 1 —p{n). 
Note that a similar model without memory was already 
defined in the literature l43l. 
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FIG. 5. Schematic definition and numerical results of the 
model study, (a) P{E) distributions of the synthetic sequence 
after logarithmic binning with window sizes At = 1...1024. 
The fitted power-law function has an exponent /3 = 3.0. 
(b) Transition probabilities of the reinforcement model with 
memory, (c) Logarithmic binned inter-event time distribution 
of the simulated process with a maximum inter-event time 
^max _ -j^Q6 rpj^^ corresponding exponent value is 7 = 1.3. 
(d) The average logarithmic binned autocorrelation function 
with a maximum lagr^"^ = 10^. The function can be charac- 
terized by an exponent a = 0.7. Simulation results averaged 
over 1000 independent realizations with parameters /xa = 0.3, 
IJiB — 5.0, V — 2.0, TT = 0.1 and T — 10^. For the calculation 
we chose the maximum inter-event time t^"^ = 10^, which is 
large enough not to influence short temporal behavior, but it 
increases the program performance considerably. 

The numerical results predicted by the model are sum- 
marized in Fig|5] and Table [l[ We find that the inter- 
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event time distribution in FigjSjc reflects strong inho- 
mogeneities as it takes the form of a scale-free function 
with an exponent value 7 = 1.3, satisfying the relation 
7 = Ma + 1- As a result of the heterogeneous tempo- 
ral behavior with memory involved, we detected spon- 
taneously evolving long temporal correlations as the au- 
tocorrelation function shows a power-law decay. Its ex- 
ponent a = 0.7 (see Figjsld) also satisfles the relation 
a + 7 = 2 (see SI). The P(E) distribution also shows fat- 
tailed behavior for each investigated window size ranging 
from At = 1 to 2^^ (see Fig|5]a). The overall signal here 
is an aggregation of correlated long bursty trains and 
uncorrelated single events. This explains the weak At 
dependence of P{E) for larger window sizes, where more 
independent events are merged with events of correlated 
bursty cascades, which induces deviation of P{E) from 
the expected scale-free behavior. The P{E) distributions 
can be characterized by an exponent /3 = 3.0 in agree- 
ment with the analytical result in ([5| and it conflrms the 
presence of correlated bursty cascades. In addition, even 
if we flx the value of /3 and 7, the a exponent satisfles the 
condition a < 7 < /3, an inequality observed in empirical 
data (see Table |l|. 



II. DISCUSSION 



In the present study we introduced a new measure, the 
number of correlated events in bursty cascades, which 
detects correlations and heterogeneity in temporal se- 
quences. It offers a better characterization of correlated 
heterogeneous signals, capturing a behavior that cannot 
be observed from the inter-event time distribution and 
the autocorrelation function. The discussed strongly het- 
erogeneous dynamics was documented in a wide range 
of systems, from human dynamics to natural phenom- 
ena. The time evolution of these systems were found to 
be driven by temporal correlations that induced scale- 
invariant distributions of the burst lengths. This scale- 
free feature holds for each studied systems and can be 
characterized by different system-dependent exponents, 
indicating a new universal property of correlated tempo- 
ral patterns emerging in complex systems. 

We found that the bursty trains can be explained in 
terms of memory effects, which can account for the het- 
erogeneous temporal behavior. In order to better under- 
stand the dynamics of temporally correlated bursty pro- 
cesses at single entity level we introduced a phenomeno- 
logical model that captures the common features of the 
investigated empirical systems and helps us understand 
the role they play during the temporal evolution of het- 
erogeneous processes. 



III. MATERIALS AND METHODS 

Data processing. To study correlated human behav- 
ior we selected three datasets containing time-stamped 
records of communication through different channels for 
a large number of individuals. For each user we extract 
the sequence of outgoing events as we are interested in 
the correlated behavior of single entities. The datasets 
we have used are as follows: (a) A mobile-call dataset 
from a European operator covering ~ 325 x 10^ million 
voice call records of ~ 6.5 x 10^ users during 120 days 
[5 . (b) Text message records from the same dataset con- 
sisting of 125.5 X 10^ events between the same number of 
users. Note that to consider only trusted social relations 
these events were executed between users who mutually 
called each other at least one time during the examined 
period. Consecutive text messages of the same user with 
waiting times smaller than 10 seconds were considered 
as a single multipart message [45 though the P(tie) and 
A{t) functions do not take values smaller than 10 sec- 
onds in Figj2]b. (c) Email communication sequences of 
2, 997 individuals including 20.2 x 10^ events during 83 
days [24]. From the email sequence the multicast emails 
(consecutive emails sent by the same user to many other 
with inter-event time 0) were removed in order to study 
temporally separated communication events of individ- 
uals. To study earthquake sequences we used a catalog 
that includes all earthquake events in Japan with mag- 
nitude larger than two between 1986 and 1998 |T2j. We 
considered each recorded earthquake as a unique event 
regardless whether it was a main-shock or an after-shock. 
For the single station measurement we collected a time 
order list of earthquakes with epicenters detected at the 
same region [3 [11] (for other event collection methods see 
SI). The resulting data consists of 198, 914 events at 238 
different regions. The utilized neuron flring sequences 
consist of 31,934 outgoing flring events of 1,052 single 
neurons which were collected with 2 millisecond resolu- 
tion from rat's hippocampal slices using fMCI techniques 
1471148.. 

Inter-event time shuffling of real sequences. For 

the independent event sequences in Figj2] and [3] (empty 
symbols) we shuffled the inter-event times of individuals 
allowing to change the inter-event time values between 
any users but keeping the original event number for each 
individual. The presented P{E) distributions were cal- 
culated with one At window size to demonstrate the ex- 
ponential behavior of P{E) for independent events. 
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IV. SCALING OF THE AUTOCORRELATION FUNCTION IN HETEROGENEOUS INDEPENDENT 

PROCESSES 

As it was mentioned in the main text for strongly inhomogeneous temporal sequences of independent events the 
Hurst exponent can take a value H > 1/2 and can assign false positive temporal correlations [1 . This effective 
behaviour is also reflected by the autocorrelation function which can show power-law scaling A{t) ^ for sequences 
of independent events with inter-event time distribution P{tie) ^ and can indicate presence of non-existing 
correlations. Based on the generating function method and the Tauberian theorems one can show that for 1 < 7 < 2 
the scaling law 

a + 7 = 2 (7) 

holds [2^. We checked this relationship by numerical simulations, where we generate independent events with inter- 
event times sampled from a power-law distribution with a fixed exponent 7 and calculate the autocorrelation function. 
As it is shown in Figjoja the exponent relation in Eqj7| holds for numerical results since for P(t^e) with 7 = 1.5 (blue 
symbols and straight line) the effective autocorrelation function (red symbols) scales as a power-law with an exponent 
a = 0.5 (dashed line). At the same time the P{E) distribution calculated with At = 10 window size (green symbols) 
indicates the true uncorrelated behaviour as it shows exponential decay. 




tie,T,E T E 

FIG. 6. The characteristic functions calculated for heterogeneous independent signals, (a) P(tie), ^(^) and P{E) 
functions for 7 = 1.5. Solid line is a power-law function with the given 7 exponent value, while dashed line denotes a a power- 
law function with an effective a = 0.5 exponent value, (b) A{t) effective autocorrelation functions for various 7 exponents. 
Straight lines are denoting power- law functions with a exponents satisfying the a + 7 = 2 relation, (c) Corresponding P{E) 
distributions for various 7 exponents. 

The same calculations had been repeated for various 1 < 7 < 2 values averaged over 1000 independent realization as 
we present it in Figj6]b. In each cases the autocorrelation function satisfies the condition derived in Eq{TT] However, 
as 7 ^ 1 extreme fluctuations start to influence the dynamics considerably, while some discrepancy appears as 7 ^ 2 
where fully correlated behaviour should evolve which cannot be the case for random processes even the fluctuations 
are finite. 

The corresponding P{E) distributions show exponential behaviour as it is demonstrated in Figjojc and assign the 
true uncorrelated temporal behaviour of the processes. It demonstrates that autocorrelation is unable to address 
present correlations obviously for heterogeneous temporal processes since it indicates effective correlations between 
independent events. However, the P{E) distribution is capable to detect correlated behaviour even for processes with 
fat-tailed inter-event time distributions as it decays exponentially in case of independent signals (for a detailed study 
see SI Section M) while it scales as a power-law when long bursty periods evolve as a result of temporal correlations. 
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V. P(E) DISTRIBUTON IN INDEPENDENT MODELS 

We study the functional behavior of the P{E) bursty event number distribution for processes of independent events. 
As it was already discussed in the main text the P{E) distribution with a given At time window size can be written 
as: 



P{E = n) 




P{tie)dU 



At 



P{tie)dU 



(8) 



forn > in case of processes where we sample Ue inter-event times independently from an arbitrary P{tie) distribution. 



If the integral P(tie)dtie < 1, then P[E = n) ^ a~^^~^^ is decreasing exponentially, otherwise if P(tie)dti, 
then P{E = n) =0. Since we fix the upper limit of the integrand At < oo, the first behavior holds. 

To numerically confirm this analytical results we define a single user model where we generate events with inter-event 
times sampled from two different distributions: 



P{Ue) - t- 



and 



P{tie) 



-t/r 



and count the number of consecutive events which fall into a bursty periods with tie ^ ^t. 



(9) 




FIG. 7. P(E) distributor in independent models. Distribution of number of bursty events in periods evolved in independent 
processes with (a) power-law and (b) exponential inter-event time distributions. Simulation results are presented with various 
parameter values and with fixed At and maximum inter-event times. 



The simulation results in Figj?] demonstrates the predicted exponential decay for the P{E) distributions. We 
performed sequences with various 7 or r values, but the P{E) distribution was calculated in each cases with a fixed 
At window size. 



VI. At TIME WINDOW DEPENDENCE OF P{E) DISTRIBUTION 
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FIG. 8. Illustration of bursty period detection in the call sequence of a selected individual. Black spikes denote the 
call events with width rational to the call length. In the 2-5 lines we colored the inter-event periods if they were smaller then 
the corresponding time window size. As we are increasing At the periods which were detected with smaller time- windows are 
be merged together to longer trains. The first line demonstrates the fine grained self-containing structure of the long evolving 
bursty periods. 
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We have seen in the main text that the P{E) bursty train size distribution shows a robust behaviour against the 
selected At time windows. By increasing At the correlated event clusters are growing as more and more single events 
and periods with shorter time windows get merged together as depicted in Fig|8] Looking for a wider range of At 
from 1 second (the time resolution of the mobile-call sequences) up to T = 120 days (the length of the data sequence) 
we can follow how P{E) distribution evolves. In Figj9]a we present P{E) distributions with At = 2^ where n goes 
from 1 to 24. A lower t[ and a higher characteristic time can be detected by looking at the evolution of P{E): 

• if At < t^: only small bursty trains evolve and the P{E) distributions show a concave behaviour. 

• if t^ < At < tjf: critical regime where the P{E) distributions are scaling as a power-law with an exponent f3 
insensitive for At. 

• if At > tj!: Uncorrelated periods are merged together and P{E) is approaching the strength distribution as 
At ^ T days where for each user only one bursty period evolves containing all events of the present user. 

For the mobile call sequence the corresponding exponent takes /3 = 4.1 and the two characteristic times are tj, = 20 
seconds and = 12 hours. These times are also play crucial roles for P{tie) and A(t) as 20 seconds is the typical 
reaction time between two consecutive call actions, while 12 hours reflects the length of correlated periods of human 
daily activity. Looking at the scaling behaviour of P{tie) and A{t) of other human activities, earthquakes or neuron 
spike sequences, one can detect analogous evolving behaviour of P{E) with corresponding characteristic times. 




FIG. 9. Distribution of correlated bursty event numbers with various time- window size. P{E) distribution calcu- 
lated for various At window size with value varying between At = 2 ^...2^^. Straight lines are denoting the lower characteristic 
time tc = 20 seconds (solid line) and the higher characteristic time tc = 12 hours (dashed line) and assign power- law functions 
with exponent (3 A.l. 



To exclude the possibility that the broad strength distribution causes the power law like P{E) behaviour, we 
repeated the same measurement with a single event sequence which was constructed from the call sequences of 10^ 
users and monitor how P{E) changes as we increase At from its minimum to its maximum (not shown here). In this 
case we also observed a extended critical regime where P{E) presented scale-free behaviour with the same exponent 



VII. STRENGTH DECOMPOSITION OF THE CHARACTERISTIC FUNCTIONS 



We demonstrated in the main text that bursty correlated behavior can be characterized by functions as P(tie), 
A{t) and P{E) in several kind of dynamic systems. However, the question remained whether the observed scaling 
behaviour of these functions is an artifact of other present inhomogeneities e.g. the broad distribution of node activity 
(strength), or it de facto reflects general behavioural characteristics of entities. To answer this question we completed 
additional measurements on the mobile-call dataset which possesses large enough size to provide good statistics after 
we decompose users into different groups. We also repeated the following measurements for the other datasets and 
found similar behaviour of the analyzed functions. 

Here we define node strength s as the number of (in and out) calls made by a user during the entire period. It 
is known from earlier studies [4 that the P{s) strength distribution of the MOD is broad with a maximum value 
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FIG. 10. Characteristic functions of different strength groups, (a) Inter-event time distributions measured between 
the outgoing calls of individuals belonging to the same strength group (see text). The solid line assigns a power-law function 
with exponent 7 = 0.7. (b) The same curves scaled by using the (t|e)^(t|e) — P (tie/ (tie)) where tfe denotes the inter-event 
times of individuals with a given s strength and (tfe) assigns the average inter-event time calculated for users with strength s. 
(c) Autocorrelation function averaged for users belonging to different strength groups. Characteristic exponents are changing 
between a = 0.7 (solid line) and a = 0.38 (dashed line) as the strength of the corresponding groups is increasing, (d) P{E) 
distribution for users belonging to different strength groups. The corresponding exponents are changing between /3 = 2.45 
(dashed line) and 7 = 4.3 (solid line). 



Smax = 7933 in the present case. Take into consideration the inhomogeneous P{s) distribution we ranked users into 
strength groups with increasing bin size. We realized from individual level analysis that users with the largest strength 
values are playing a disparate role as they show non- human like communication patterns. Though we exclude them 
from the measurements and take users only with strength smaller then 96% of the maximum value. 

In Fig{T0]we present the average characteristic functions calculated separately for each strength groups. Hetero- 
geneous inter-event time distributions are characterizing the communication of users ranked into different groups 
(Fig|lQ|a). Only the tail of the P(t^e) distributions show significant discrepancy due to different activity levels. Since 
users frequently show correlated bursts with short inter-event times, those ones with small number of calls in sum 
have longer inactive periods between bursts which induces a longer tail in P{tie) with a later exponential cutoff. Users 
with higher activity have shorter waiting times between bursts which is reflected by the earlier turning point of the 
inter-event time distribution. In order to proof whether these distributions are broad not only due to the inhomoge- 



neous strength distributions in Fig 10 b we scaled them together using the scaling relation (^|e)^(^|e) = ^(^le/(^le)) 
where (tf^) denoted the average inter-event times, calculated separately for each groups as it was done in [5, 6 . Using 
this scaling relation the distributions scale together on the same master curve which indicates that P{tie) follows the 
same functional behaviour independently from the chosen strength group (and the average inter-event time of this 
group). 

The autocorrelation function in FigfTolc decays as a power-law up to 12 hours assigning long-temporal correlations 
for users in each strength groups. Naturally stronger correlations are detected for more active users as the corre- 
sponding exponents vary between a 0.38... 0.7 as we decrease the activity level. A similar behaviour is confirmed 



by the scaling of P{E) in Fig 10 d as the distributions (calculated with At = 600 seconds) remain fat-tailed for each 
user groups with an exponent between P = 2. 45... 4. 3 as we decrease activity. It implies that long correlated bursty 
periods evolve even for users with only a few calls but with smaller probability then for users with many call actions. 
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VIII. DE-SEASONED RESULTS 



In order to study the effect of circadian patterns on the observed behaviour of the characteristic functions in mobile 
cah communication, we remove daily fluctuations by rescaling the event times using a method defined in [71 E]. 
Measuring the r{t) event density and its average value Rt during the entire period we can define a rescaled time for 
each event as: 

(it* = ^dt = p{t)dt (10) 

where p{t) denotes the event rate and the p''{t*)dt* = p{t)dt scaling holds for the time variable with p''{t*) = 1. 
Rescaling event times with these conditions, events in periods with high event frequency become dilated, while event 
times are contracted when their frequency is lower, so the effect of cyclic fluctuations can be reduced. We rescaled 
the times of the outgoing call events of each user in the MCD considering the duration of their calls and ranked them 
into strength groups using their overall call activity. 




FIG. 11. Characteristic functions of de-seasoned sequences, (a) Autocorrelation functions and (b) inter-event time 
distributions of de-seasoned outgoing call event sequences of individuals belonging to the same strength groups. Dashed line on 
panel (a) assigns a power-law function with exponent a = 0.65 while on panel (b) with exponent 7 = 0.7. (c) P{E) distributions 
measured with various At time window size in the de-seasoned outgoing event sequences of individuals. The slope of the dashed 
straight line indicates an exponent value /3 = 4.1. (d) P{E) distribution with At = 600 measured for users belonging to different 
strength groups. The solid (dashed) line denotes a power-law function with exponent 4.6 (2.8.) 



Utilizing the sequence of outgoing call events with rescaled times we calculate the three characteristic functions to 
see the impact of daily fluctuations on them. Since cyclic fluctuations are introducing correlations in the dynamics, by 
removing them the A{t) function should reflect weaker correlations and if only the circadian patterns are responsible 
for the temporal correlations this function should radically change and present short term correlated behaviour only. 
In FigfTHa the autocorrelation denotes reduced correlated behaviour compare to the same function of the original 
event sequence in Fig 10 c, however it signifies remaining long temporal correlations as it shows slow decay with a 
slightly larger exponent a o:^ 0.75 compared to the unsealed value {a = 0.5). The inter-event time distributions in 
Figjnjb also remains similar compared to FigjlOja with approximately the same exponent value 7 0.7. It is in 
complete agreement with the results presented in [7] where the P(tie) of call events remained unchanged after the 
same de-seasoning method was applied on the event sequence. 

Long bursty periods are not destroyed by removing daily patterns from the event sequence as it is shown in 
Fig(TT|c. The long periods of outgoing bursts of individuals remain for various At values reflected by power-law like 
P{E) distributions decreasing with exponent ^ A.l similar to the original sequence. The P{E) distribution becomes 
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more disperse if we decompose it by strength for periods with time window size At = 600 seconds (see FigjTTjd). 
The tail of the decomposed P{E) distributions can be estimated with exponents between /3 2. 8... 4. 6. Consequently 
as the scaling behaviour of the de-seasoned data show similar behaviour as the original sequence, it implies that 
even circadian fluctuations are partially responsible for the present correlations they do not effect considerably the 
evolution of long correlated bursty sequence in communication of individuals. 



IX. BURSTY-TOPOLOGY CORRELATIONS IN EARTHQUAKE-SEQUENCES 



In order to measure temporal correlations between earthquake events we applied a so-called single-station method 
[OHTlj and studied event sequences executed at the same geographical area. For each earthquake event the longitude 
and latitude coordinates of the epicenter was known from informations measured at several surrounding seismological 
stations fT2|. The distribution of epicenters on the investigated area is visualized on Fig 12 a. 
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FIG. 12. (a) Earthquake events in Japan between 1985.07-1998.12. Count of earthquake events located into area bins with size 
(b) 10°A:m^ (c) lO^/cm^ and (d) lO^km^. 



For the results presented in the main text we grouped 198, 914 events into 238 administrative regions and we observed 
scale-free behaviour for the characteristic functions (see Fig. 3. a in the main text). However, since the regions can have 
different sizes the question remains whether the observed scaling behaviour is the result of the diverse size distribution 
of different regions or it truly assigns heterogeneous correlated temporal behaviour. To answer this question we divided 
the investigated area for bins with equal sizes and for each bin we collected a time-ordered list of events executed on 
the corresponding area. We used three different bin sizes: 10^ /cm^, 10^ /cm^ and 10^/cm^. The event count map with 



the three different bin sizes is shown in Fig 12 b, c and d. 

Using the event sequences collected for each bins we calculated the characteristic function P(t^e), ^(t) and P{E). 
As it is shown in Fig 13 each function remained fat tailed with the same exponent values P = 2.5 and 7 = 0.7 as 
it was found for the measurements presented in the main text. The autocorrelation function exponent took a value 
a = 0.25 which assigns slightly stronger correlations as for the earlier calculations, however the average autocorrelation 
functions here were calculated for the first 100 most active bins which explains the slight decrease of the a exponent. 



10-^ ► 
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FIG. 13. Characteristic functions: (a) P{tie), (b) A{r) and (c) P{E) calculated for sequences collected with area bin sizes 
10° /cm^ (blue), lO^/cm^ (red) and lO^km^ (grey). The autocorrelation functions on (b) were averaged for the first 100 most 
active bins with the corresponding size. On (c) the different symbols denote the P{E) distributions calculated with different 
At time window sizes. 



As the scaling behaviour of the characteristic functions were re-found for calculations with equal-size area bins, in 
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conclusion we can say that the observed temporal correlations and heterogeneous dynamics are not the result of the 
event collection method but they truly characterize the sequences of earthquake events. 



X. CALCULATION OF THE MEMORY FUNCTION 



In the main text for correlated event sequences with P{E) ^ E ^ we can derive the memory function in the form: 



Pin) ={—,] (11) 



However, by only knowing EqjTTJwe can also calculate the related P{E) distribution as: 

---('-(;^)Tn:(-,)^(.-(-,r)a)^ 

_ (n + If -n^ + n-^y - 1 _ + n + + -j " 1 



n{n+lY n^'+i 



(13) 
(14) 



Consequently our calculations are consistent as the corresponding P{E) distribution shows asymptotically a power-law 
decay with an exponent satisfying the relation: 

/3 = z/ + l. (15) 



XL ESTIMATING EXPONENT FOR EMPIRICAL MEMORY FUNCTIONS 




FIG. 14. Fitted complement memory functions of different empirical data sequences calculated with a selected 
time-window size, (a) Mobile call sequence (At = 600 seconds), (b) Text message sequence (At = 300 seconds), (c) Email 
sequence (At = 3600 seconds), (d) Japanese earthquake sequence (At = 16 hours), (e) Neuron firing sequence (At = 900 
milliseconds) . 



In order to estimate the v memory function exponents for empirical sequences we calculated the p{n) function for 
each investigated data sequence with a chosen At time window size and plot the complement 1 — p{n) of the original 



16 



memory functions in Fig {m] (grey symbols) and the same functions after logarithmic binning (black symbols). As it 
is visible the original complement memory functions show strong finite size effects as the investigated event sequences 
span on a limited time frame and the related P{E) distributions are finite. Similarly as it was discussed in the 
main text, we fitted the binned empirical memory functions with the analytical functions of the form of EgfTT] using 
least-squares method with only one free parameter, the exponent u (red lines in Fig 14). The resulted u exponents are 
written in the figure caption and also summarized in the main text in Table I. The derived exponents approximately 
satisfy the relation in Eq{T5]with the corresponding empirical /3 values. 

The actual p{n) memory functions take the analytical form Eq TT]if we assume that the P{E) distribution is a power- 
law function with an exponent 7. However, the fitted curves show deviations from the empirical p{n) functions for 
small n values as the measured memory functions are derived from non-perfect empirical power-law P{E) distributions. 
Nevertheless, the analytical and empirical functions fit well asymptotically for larger n values. 
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